import ps_util as ps
import numpy as np
import matplotlib.pyplot as plt
import plotly
plotly.offline.init_notebook_mode()
%matplotlib inline
#read in the file which has I,M,S where I is the images, M is the mask, and S is the light vectors.
V=ps.read_data_file('Beethoven')
# Take the I,M,S and store them into variables
Iv=V[0]
mv= V[1]
sv= V[2]
plt.figure(figsize=(7,7))
plt.imshow(Iv)
plt.show()
plt.figure(figsize=(7,7))
plt.imshow(mv)
plt.show()
plt.figure(figsize=(7,7))
plt.imshow(sv)
plt.show()
#displaying the images
for i in range (0,3):
plt.figure(figsize=(7,7))
plt.imshow(Iv[:,:,i])
plt.show()
# use the mask_image to calculate the relevant values within the mask (where values > 0)
nz=np.where(mv)
len(nz[0])
28898
mv.shape
(256, 256)
# Take each individual image and apply the mask to it, then reshape the outupt into a (3,nz) array.
a=Iv[:,:,0]
b=Iv[:,:,1]
c=Iv[:,:,2]
ai=a[nz]
bi=b[nz]
ci=c[nz]
J=np.row_stack((ai,bi,ci))
J.shape
(3, 28898)

I is the intensity of each pixel in the image, p is the albedo of each pixel, N is the surface normal and S is the lighting source with the assumption that it is constant over each image.
Take the dot product of the Image times the inverse of the light vector to obtain the Modulated Normal Field.
#take the light intesnity vector, take its inverse and store it in a variable to use in Lamberts Equation
# M = S^-1 * J
SV=np.linalg.inv(sv)
SV.shape
print(SV)
[[ 1.62948467 -3.51636848 1.70279567] [ 1.95828261 -0.13769666 -1.75540327] [ 0.36274532 0.32154769 0.31388367]]
#Take the dot product of the two matrices in order to solve for M
Mb=np.dot(SV,J)
Mb.shape
(3, 28898)
Albedo: material light absorption property
#Normalize the M to obtain the albedo
albedo=np.linalg.norm(Mb, axis=0)
albedo.shape
(28898,)
#Store the albedo back into an empty numpy array and visualize the output.
Beto=np.zeros((256,256))
Beto[nz]=albedo
fig, ax = plt.subplots(figsize=(7,7))
plt.imshow(Beto)
plt.title("Albedo for Beethoven in 2D")
plt.axis('off')
(-0.5, 255.5, 255.5, -0.5)
#Albedo with a different color scheme
plt.figure(figsize=(7,7))
plt.imshow(Beto, cmap='binary')
<matplotlib.image.AxesImage at 0x10ef03a90>
#Calculate the normal fields n1, n2, n3
n1=np.divide(Mb[0],albedo)
n2=np.divide(Mb[1],albedo)
n3=np.divide(Mb[2],albedo)
Nbus=np.vstack((n1,n2,n3))
#Take the normals and store them into an empty numpy
be1=np.zeros((256,256))
be2=np.zeros((256,256))
be3=np.ones((256,256))
be1[nz]=n1
be2[nz]=n2
be3[nz]=n3
Betov=np.dstack((be1,be2,be3))
Betov.shape
(256, 256, 3)
#code to use to plot the image through plotly
x = np.where(Betov)[0]
y = np.where(Betov)[1]
x.shape
y.shape
(123332,)
#display images
plt.figure(figsize=(6,6))
plt.imshow(be1)
plt.title("Normal Field 1")
plt.show()
plt.figure(figsize=(6,6))
plt.imshow(be2)
plt.title("Normal Field 2")
plt.show()
plt.figure(figsize=(6,6))
plt.imshow(be3)
plt.title("Normal Field 3")
plt.show()
plt.figure(figsize=(6,6))
plt.imshow(Betov)
plt.title("All Normals")
plt.show()
#display images
plt.figure(figsize=(6,6))
plt.imshow(be1, cmap='binary')
plt.title("Normal Field 1")
plt.show()
plt.figure(figsize=(6,6))
plt.imshow(be2, cmap= 'binary')
plt.title("Normal Field 2")
plt.show()
plt.figure(figsize=(6,6))
plt.imshow(be3, cmap='binary')
plt.title("Normal Field 3")
plt.show()
plt.figure(figsize=(6,6))
plt.imshow(Betov, cmap='binary')
plt.title("All Normals")
plt.show()
#calculate the z to plot using unbiased integrate
z0=ps.unbiased_integrate(be1,be2,be3,mv)
z0.shape
(256, 256)
ps.display_depth_matplotlib(z0)
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
#calculate the z to plot using simchony
z1=ps.simchony_integrate(be1,be2,be3,mv)
ps.display_depth_matplotlib(z1)
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:584: RuntimeWarning: invalid value encountered in less
cbook._putmask(xa, xa < 0.0, -1)
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:1545: RuntimeWarning: invalid value encountered in greater
intensity > 0),
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:1550: RuntimeWarning: invalid value encountered in greater
hsv[:, :, 2] = np.where(intensity > 0,
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:1556: RuntimeWarning: invalid value encountered in less
intensity < 0),
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/colors.py:1560: RuntimeWarning: invalid value encountered in less
hsv[:, :, 2] = np.where(intensity < 0,
/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
if self._edgecolors == str('face'):
#since matplotlib and mayavi were subpar options, we opted to graph using a plotly interactive graph
import plotly.plotly as py
from plotly.graph_objs import *
trace1 = Surface(
z=z1, # link the fxy 2d numpy array
x=y, # link 1d numpy array of x coords
y=y, # link 1d numpy array of y coords
colorscale="Greys"
)
data = Data([trace1])
# Dictionary of style options for all axes
axis = dict(
showbackground=True, # (!) show axis background
backgroundcolor="rgb(204, 204, 204)", # set background color to grey
gridcolor="rgb(255, 255, 255)", # set grid line color
zerolinecolor="rgb(255, 255, 255)", # set zero grid line color
)
# Make a layout object
layout = Layout(
title='Beethoven', # set plot title
scene=Scene( # (!) axes are part of a 'scene' in 3d plots
xaxis=XAxis(axis, showticklabels=False), # set x-axis style
yaxis=YAxis(axis, showticklabels=False), # set y-axis style
zaxis=ZAxis(axis, showticklabels=False) # set z-axis style
)
)
# Make a figure object
fig = Figure(data=data, layout=layout)
#the graph is interactive and you can zoom in and out of the 3D figure
plotly.offline.iplot(fig)